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Calculation of moments of generalized parton distributions in lattice QCD requires more powerful 
techniques than those previously used to calculate moments of structure functions. Hence, we present 
a novel approach that exploits the full information content from a given lattice configuration by 
measuring an overdetermined set of lattice observables to provide maximal statistical constraints on 
the generalized form factors at a given virtuality, t. In an exploratory investigation using unquenched 
QCD configurations at intermediate sea quark masses, we demonstrate that our new technique is 
superior to conventional methods and leads to reliable numerical signals for the n = 2 fiavor singlet 
generalized form factors up to 3 GeV^. The contribution from connected diagrams in the fiavor 
singlet sector to the total quark angular momentum is measured to an accuracy of the order of one 
percent. 
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I. INTRODUCTION 



Light-cone correlation functions play a special role in the experimental exploration of the quark and gluon structure 
of hadrons. Asymptotic freedom allows quantitative separation of the reaction mechanism from the structure of 
the probed hadron at high energy, so that spin-independent scattering experiments unambiguously measure matrix 
elements of the light-cone operator 

0{x) = / ^e^^^^(-^n) /i7'e-'^-f-V2^«"-^("")^(^n), (1) 
J in 2 2 

where n is a light-cone vector, and V denotes a path-ordering of the gauge fields in the exponential. Since these matrix 
elements are singled out by their experimental accessibility, it is essential to use all our tools of analytical methods 
and lattice field theory to explore and understand them as fully as possible. 

Diagonal nucleon matrix elements, q{x) = {P\0{x)\P), measure the familiar quark distribution q{x) specifying the 
probability of finding a quark carrying a fraction x of the nucleon's momentum in the light cone frame. Although 
light cone correlation functions cannot be calculated directly in lattice QCD, expansion of 0{x) generates the tower 
of twist-two operators, 

0^^i'^=-''"> = ^^-f^^''iD^"- . . . lD^'"^g , (2) 

with = ^ (Id and curly braces {/ii . . ./i„} mean symmetrization of indices and subtraction of traces. The 

diagonal matrix elements (pjo^'^i^^.-.A'.} 

can be calculated on the lattice and specify the {n — 1)*'* moments 
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/ dx x'^~^q{x). Note that expressions analogous to Eqns. ^ and (|2Jl for spin-dependent observables differ only in 
their gamma matrix structure, but will not be considered in the present work. 

Generalized parton distributions (GPDs), as introduced in correspond to non-diagonal matrix elements 

{P'\0{x)\P) . When expressed in terms of the relevant Lorentz invariants, {P'\0{x)\P) is specified in terms of two 
generalized parton distributions, H{x,^,t) and E{x,^,t), depending on three kinematical variables. In terms of 
the four- momentum transfer A = P' — P, the invariant momentum transfer squared is t = A^, the skewedness is 
^ — —n ■ A/2, and x denotes the momentum fraction. Since the dependence of the GPDs, H{x,£,,t) and E{x,£,,t), 
on three kinematical variables renders their physical interpretation more difficult than ordinary parton distributions, 
it is useful to recall several important physical properties. In the forward limit, i.e. ^, t — > 0, we recover the 
forward parton distribution function as H[x,0,0) = q{x). In what is sometimes referred to as the local limit, 
integrating over the momentum fraction, x, yields the familiar electromagnetic form factors, J dx H{x,^,t) = Fi{t) 
and J dx E(x,S^,t) — F2{t)- The first moment of the sum of H and E yields the total quark angular momentum 

J dx x[H{x, ^, t) + E{x, t)] — Jq. We note that both these results are independent of ^. And finally, at skewedness 
f = 0, the t dependence specifies the transverse Fourier transform with respect to the impact parameter of the light 
cone wave function j^]. For a possible interpretation of GPDs with non-zero ^, see Refs. 5, ,6j. 

On the lattice, instead of matrix elements of the light cone operator, Eq. one again calculates non-diagonal 
matrix elements of the local operators, Eq. (|2Jl, yielding moments of the generalized parton distributions. Following 
the notation of Ref. the non-diagonal matrix element (P'|C'f^i- -^'"i|P) rnay be expressed in terms of a set of 
generahzed form factors (GFFs) Ani{t), Bni{t), and Cn{t). The form factors Ani{t) and C„(t) multiphed by powers 
of yield the moment Hn{£,,t) = J dx x'^^^H{x,£,,t) and the form factors Bni{t) and Cn{t) multiplied by powers of 

yield the moment En{S,,t) — J dx x'^~^E{x, ^,t). 

The lowest three moments considered in this work are 

{p'\o^^\p) = {(riA.oit) 

+ ^(K^"»A„i?ioW, 

+ -^Ft''i((a''^>"))A„i?2oW 
2m 

+lA{MiAM2}((i))C2(i), 
m 

+A^^iA^^((7^^>))A32(0 
2m 

where = (P^ + P'^)/2 and ((F)) = U{P')TU{P). The GFFs, P„,(i), and C„(t), specify all the information 

about spin-independent generalized parton distributions that is known to be accessible on the lattice. The limits for 
H and E discussed above may be re-expressed in terms of the generalized form factors. The limit t ^ {) oi A„o is 
the familiar parton distribution moment, yl„o(0) = J dx x"'~^q{x). The electromagnetic form factors are given by 
^10 (t) — Fi {t) a-iid BiQ (t) = F2 (t) for the appropriate flavor combination. Finally, the total quark angular momentum 
is given by the sum of A20 + P20 as i ^ 0, Jg = ^[^20(0) + -820(0)]. 

In the context of this brief review, we may now consider the challenges and opportunities in calculating moments 
of generalized parton distributions on the lattice, and compare them with the analogous issues for ordinary parton 
distributions. Experimentally, three decades of deep inelastic scattering experiments haveprovided impressive phe- 
nomenological determinations of parton distributions as a function of momentum fraction j^ilSllSl- Hence, for parton 
distributions, the key issues are developing the lattice technology to the point of attaining quantitative agreement with 
moments of experimental parton distributions and using the lattice as a tool to obtain insight into how these distri- 
butions arise from QCD. The present status is that computational limitations restrict unquenched QCD calculations 
to the heavy quark regime in which the pion mass is heavier than roughly 500 MeV and naive linear extrapolation 
to the physical pion mass yields serious disagreement with experiment. For example, linear extrapolation of the 
quark momentum fraction exceeds experiment by the order of fifty percent 0, 0, 0, 0, ■ Although there are 
strong indications that physical extrapolation to the chiral limit may introduce corrections of the required magnitude 
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|15|. there is presently no quantitative theory for the intermediate mass regime and we must exploit the emerging 
generation of computers to perform the requisite calculations sufficiently close to the chiral regime. 

For generalized parton distributions, the situation is quite different. Since GPDs depend on three variables and 
experimental quantities involve convolutions, there is no prospect of measuring the full dependence on x, ^, and 
t. Without additional input arising from first principles, extraction of GPDs from experiments like deeply virtual 
Conipton scattering will necessarily be contaminated by uncontrolled assumptions. Hence, once computer power is 
sufficient to obtain quantitative agreement with moments of parton distributions, lattice calculations of moments of 
GPDs will become an essential tool to be used in conjunction with experiment to extract and understand the full 
dependence on cc, ^, and t. It is thus imperative to develop techniques to calculate these moments. 

In addition, theorists may also obtain insight into how QCD works by studying the dependence of hadron structure 
on the quark mass. This study can begin immediately, addressing the behavior of hadrons in a world where the pion 
weighs more than 500 MeV. This heavy pion world is much closer to the non-relativistic quark model, and as we 
eventually lower the pion mass, we will learn how QCD evolves from the world of heavy quarks to the physical world 
of light quarks. Even in the heavy pion world, one can test contemporary assumptions, such as factorization of the t 
dependence [l^ . 

On the computational side, it is essential to confront the additional challenges that arise for GPDs relative to the 
ordinary parton distributions. Already for forward parton distributions, the tower of operators, Eq. (0), involves 
operators that become increasingly subject to statistical noise as one progresses to higher and higher derivatives. 
For GPDs, however, we compound the noise of these operators with the additional noise from the finite momentum 
transfer, A. Hence, in this work, we address the problem of imposing the maximal statistical constraints a lattice 
calculation can provide on the form factors Ani, Bni, and C„ appearing in Eq. lO. We note that we have at our 
disposal the choice of several alternative representations of the hypercubic group on the lattice corresponding to the 
same continuum operator 0{A'iM2. -Aire} g^g g^g choice of different kinematic variables corresponding to the same 
t. We will therefore use this freedom to construct an over-determined set of lattice observables corresponding to the 
continuum expressions in Eq. ((JJ, and thereby significantly improve the measurement of the form factors. In section 
II, we will describe the details of the method. Section III will compare the results of our method with a conventional 
analysis and present results for the the angular momentum carried by quarks, Jq. The conclusions will be presented 
in the final section and an Appendix presents the necessary detailed expressions for matrix elements in terms of 
generalized form factors. 



II. LATTICE CALCULATION OF GENERALIZED FORM FACTORS 

We will calculate the matrix elements of ^p'|0^'^i^^- -^"^ to extract the generalized form factors ^^^(t), B'^^{t), 
and C^(t) where, when relevant, we append a quark flavor label q. Since the calculation of disconnected diagrams 
raises yet another level of complexity, in this present work we will restrict our attention primarily to the flavor 
non-singlet combination u — d, for which these diagrams cancel. 

In the usual way, we calculate the non-diagonal matrix elements by the following ratio of three- and two-point 
fimctions, 

C^P^T-^nk -T + T,,e, P) C^^^T, P') C^P^T^nk, P') ] ^'^ , . 

C2pt(T,„k - r -f T,,e, P') C2pt(r, P) C2pt(r,„k, P)\ ' ^ ' 

The factors relating this ratio to the physical continuum matrix element are given in Eq. (|A2(I . The correlation 
functions are given by 

C2P*(t,P) = ^(ru„pol),,(f^|iVfc(T,P)F,(T,,e,P)|l}) , 

C^^\t, P', P) (rp°i)j-fc ^^(^-k, P')0{t)N,{t,,,,P) \n) , (5) 

where \Vl) denotes the QCD vacuum state. The nucleon source, N{t,P), and sink, N{t,P), create and annihilate 
states with the quantum numbers of the nucleon, and to maximize the overlap with the ground state, we used the 
smeared sources defined in Ref. 0|. The source is located at timeslice Tgrc, the operator O is inserted at timeslice r, 
and the sink is positioned at timeslice Tgnk- Explicit expressions for the polarized and unpolarized projectors Fpoi/unpoi 
are given in Appendix 1X1 Eqs. (|A3I A4|) . 



, , , C?^\t,P',P) 
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Inserting a complete set of states into Eq. ((SJ and using the time evolution operator yields 
C^P* (r, P) = ^ e-^' (^^) Tr (r^npoi (mir, P) I /) {I \N{nrc ,Pm) 



I 

+higher states , (6) 



- „--Eo(P)(r^r,„)--Eo(P')(r,„i.-r) {Z{P)Z{P')) JT(P')U(P\\ 

~ ' Eo{P')Eo{P) ^nh^oiU[r )U[r)) 

X (P'l 0(r) |P) + higher states . (7) 

The contributions from higher states in Eqns. lO and ((TJ are suppressed by exponential prefactors when rs„k ~ t and 
T — Tsrc are significantly greater than the inverse of the excitation energy of the first excited state. 

The ratio, Eq. Q), is constructed to exactly cancel all exponential and wave- function overlap factors. The two-point 
functions C^P'(rsnk7 P') and C^P'(Tsnk, P) decay exponentially for the full Euclidean distance between the source and 
sink and are thus particularly subject to statistical noise with finite statistics. In the worst case, they may even 
become negative, and these cases are excluded from the present work. We note that other possibilities besides Eq. 4 
may be used to cancel the exponential and overlap factors, and this freedom will be explored in a subsequent work. 

For sufficiently large time separations the ratio P(r, P', P) will exhibit a plateau yielding the desired lattice matrix 
element, and the plateau value, R{P',P), is obtained by averaging over an appropriate range of timeslices, Tmin to 

R{,^f.....f.^}{P\P)^- E Ii{f^i,2...,,M^:P':P)- (8) 

'max 'mill 

To convert our lattice calculations to the continuum MS-scheme, we use 1-loop perturbative matching at the scale 
^2 = 4GeV2 

j 

so that the lattice matrix element is related to the continuum matrix element by 

(P'lOpiP) = VEiP')EiP)Y,Z,,R,. (10) 

j 

Note that the renormalization constant, Za{fj,, a), depends only on the operator O, but not on the external states, so 
that Eq. H10|l is valid for any external momenta, P' and P. 

Finally, we write the Euclidean continuum relation between the renormalized matrix element of the generalized cur- 
rent {P'lO"^^^^^ ^ ||P) and the desired generalized form factors A'^^^{t), B'^^it), and C^(t) in the following abbreviated 
notation: 

(^'1 <^U,....,.} \p) - E + E + ■ (11) 

i 3 

Full expressions for the kinematic factors {a^, bj, c} are given in Appendix 1X1 for n ~ 1,2, and 3, Eq. ljA5|l . 

For a given n, we may evaluate the N^^^ generalized form factors {A'!l^i{t) ^ B'^-{t) , C1{t)} as follows. We select 
N > N^^^ sets of operators C^^^^^ ^ j and momenta {P', P} such that Eq. specifies N^^^ linearly independent 
combinations of the form factors. Lattice matrix elements for these operators and momenta are calculated from the 
the ratios R{t,P' ,P), in Eq. Q and matched to contimmm operators via Eq. If = N^^^, the GFFs are 
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calculated by inverting Eq. Hll|l . and if iV > N^^^ , they are calculated by a least-squares fit to the overdetermined 
system Eq. ([TT|) . 

Note, that in contrast to the case of forward parton distributions, where the moments correspond to a single number 
(denoted w„ or (a;"~^)), or electromagnetic form factors, where there are two form factors {Fi and F2), we have the 
complication of N^^^ unknown generalized form factors to be determined. We therefore now discuss the strategy for 
selecting an appropriate set of operators and momenta for this task. 

A. Practical considerations 

One practical concern is the numerical noise associated with momentum projection. The three-point function is 
subject to noise from the projection of the sink onto momentum P' and of the operator onto momentum transfer A. 
In addition to being subject to the sink momentum projection, the two-point functions appearing in the lattice ratio 
also collect noise from the projection onto the source momentum P. Ideally, for each invariant momentum transfer 
t, one would like to select P and P' such as to minimize these errors, but each distinct choice of sink momentum 
P' requires an expensive calculation of a new set of backward propagators, whereas changing A for fixed P' requires 
no new propagators. Hence, denoting the lowest momentum attained on a lattice with lattice points in a spatial 
dimension as pi = 2t:/ Ls- As a practical matter, we have used two sink momenta, P' = (0, 0, 0) and P' = {—pi, 0, 0). 
With these sinks, we generate a substantial set of momenta listed in Appendix IXI Tab. ITU by using a range of values 
of A. However there are two special virtualities of note, the virtuality, iviow and the Breit frame virtuality tercit given 
by the following momentum combinations, or rotations thereof: 

P' = (0, 0, 0) , P = (p,, 0, 0) ^ tvlow =(m- ^pf + to2^ -pf, (12) 

P'- (-pi,0,0),P=-P' ^ tBrcit=4pf, (13) 

The virtuality iviow is the lowest virtuality that can be placed on the lattice for the set of external momenta we use, 
for which no spatial momentum exceeds pi, and would be expected to have minimal projection error. Since the time 
component does not reduce the virtuality in the Breit frame, the Breit frame provides the optimal means of providing 
a large virtuality with minimal momentum projection error. Hence, including both iviow and ^Brcit in our full set 
of momentum selections ensures the presence of measurements at both ends of our t-range that have the minimal 
possible momentum-induced statistical error. In addition to these considerations, for certain different virtualities 
there are more available momentum combinations and thus more constraints than in the Breit frame. Hence, we 
obtain comparably small errors for most of the other virtualities as well. 

Another practical issue in measuring generalized form factors is the presence of powers of the momentum transfer 
A in the kinematical factors in Eq. Ijllfl . The four- momentum A is — compared to the time component of external 
momenta in the cases we consider — a small number, so its presence in kinematical factors amplifies the effect of 
statistical errors in lattice matrix elements on measurements of the associated form factors. Consistent with the 
pattern observed in Eq. |31 in the general case the factors multiplying AJ^j(i) contain i = {0,2,- • •2[^^Yi]} factors 
of A, the factors hj multiplying Bl^-it) contain j = {1, 3, • • • 1 -I- 2[^^y^]} factors of A, and the factor c, which is non- 
vanishing only for even n, has n factors of A. Therefore, we expect ^^o(0 t° quantity that can be extracted 
most accurately from lattice calculations, P,J^o(0 should have slightly larger errors, C^(t) is the worst determined when 
it is non- vanishing, and the remaining GFFs, A'^^^{t) and B'^-{t), should lie somewhere in between these extremes. 

B. Overdetermined set of lattice observables 

To extract generalized form factors from lattice calculations, we now consider a fixed value of the virtuality, that 
can be achieved on the lattice, and abbreviate Eqns. (|10() and (|11() in the following schematic form: 

{or') = E«*^-^^' 

{or') - v:b^e^.,p,, (14) 

3 

where we denote the generalized form factors generically by J^j(t) and j runs over all the form factor labels for the 
n under consideration. Thus, for n = 1, J-i{t) = AlQ{t) and J-2{t) = A2Q{t). Note that the familiar spin non-flip 
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and spin-flip electromagnetic form factors are defined by these flavor form factors weighted with the electric charge, 
i.e. Fi{t) = 2/3A^o(i) - l/3^io(i)> and F2{t) = 2/3B5^o(0 - l/3S5'o(0- Similarly, for n = 2, corresponding to the 
operator the in Eq. ^ denote J^i(i) = ^^^(i), J^2(0 = S|o(0> and J^3(t) = C|(t). 

Eliminating (O™"*) from Eqs. (|14|l yields the final form of the set of equations we will use for our analysis: 

3 

In this set of equations, the label j runs over the N^^^ generalized form factors for the given n. The label i denotes a 
specific choice of nucleon momenta P and P', and an index combination for a lattice H(4) representation of a continuum 
operator O^^^^^ ^ j- We consider N distinct values of this label, i, corresponding to the same n and t and require 

N > N^^^ . Operationally, to obtain this collection of N equations, we list all external lattice momenta, group them 
together in classes of identical virtual momentum transfer, and write down for each momentum combination the set of 
indices, {i}, of 0{i} that gives rise to non-trivial kinematic factors, aij. Note that sets of external momenta in which 
the spatial components differ only by cubic rotations still yield different constraints. Rotations changing components 
along the spin polarization axis in general change the kinematic factors and yield independent equations. Rotations 
that change components orthogonal to the spin polarization axis yield constraints that are physically equivalent, 
but provide additional statistical constraints in a lattice Monte Carlo calculation. Appendix El shows the explicit 
momentum combinations and the basis of diagonal index combinations we have chosen in this work. 

Equation (|15() is an overdetermined system of linear equations that in principle provides the maximal information 
on the generalized form factors that is attainable from the lattice. We will compare it in Sec. IllI Al with the conven- 
tional technique of using a uniq uely determined system of equations. In the case of measuring moments of parton 
distributions (see, e.g. Refs. [lOl IllL fO. Il3l Il4| l the conventional technique corresponds to selecting a single row of 
Eq. 1)15(1 . and for electromagnetic form factors, it corresponds to taking linear combinations or selecting kinematic 
variables such that Fi and F2 appear separately in two equations. 

It is important to note that Eq. (|15|l contains systematic, as well as statistical errors. The lattice action, the discrete 
approximations to the operators C'^^^^ n }' ^'^'^ ^^"^ nucleon dispersion relation contain lattice artifacts, and since 
we use Wilson fermions, the calculation contains errors of 0{a). In principle these lattice artifacts can by addressed 
by improved operators and extrapolation to the continuum limit, but these corrections are beyond the scope of the 
present work. In addition, the use of 1-loop perturbative renormalization introduces systematic errors in the analysis, 
which ultimately can be improved by higher order corrections or non-perturbative renormalization. 



C. Singular Value Decomposition 

We now turn to the solution of the overdetermined set of equations, Eq. ((15|l . which we write in matrix form as 
follows: 

R' = A'-T. (16) 

Here, T is the vector of the desired GFFs, R' denotes the statistically measured lattice ratios, and A' is the matrix of 
the coefficients a'^. Changing to a standard notation, we let n -I- 1 denote the number of GFFs, which we previously 
called N"^^^ , and continue to call the number of constraints N, so that A is an x (n -|- 1) matrix, ^ is a vector of 
length n + 1, and R' a vector of length N. 

We solve this overdetermined problem by minimization of the norm 



,2 



^.^^/E^iik:^] , (17) 



where ai denotes the jackknife error of the lattice measurement of i?^. In this case the jackknife procedure should be 
performed with an appropriate bin-size to eliminate autocorrelation effects. It is convenient to absorb the errors ai 
by defining Aij = A'^j/ai and Ri = R'i/(Ji so that we seek a vector J- that minimizes 



X^ = \A-T-R\\ 



(18) 
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The method of choice to solve Eq. H18|l in the presence of singularities or near singularities is Singular Value Decom- 
position (SVD) [13. It is based on the theorem that any x (n + 1) matrix may be decomposed as 

A = U ■ diag(u;i, . . . ,w„+i) • V , 

where J7 is an x (n + 1) matrix, V is {n + 1) x (n + I), and all Wi, the singular values of A, are non-negative. If all 
the Wi are non-zero, the solution that minimizes Eq. p8(l is given by 

F^V- diag [l/wi, l/wn+i) ■ {U^ ■ R) , (19) 

which provides the optimal vector of form factors, we can achieve with all the information available. 

If none of the singular values is zero, the method is equivalent to standard least-squares minimization. If one or more 
of the singular values, Wi, are zero, the rank of the matrix A is reduced by the number of zero singular values, implying 
that the associated directions in the solution space can not be explored. In the case of zero or nearly zero singular 
values, we may opt to avoid exploring the offending direction(s) and minimize Eq. (|18() in the residual subspace by 
replacing the corresponding factors 1/wi by zero. In this case, of course, only the components of J- corresponding to 
non-zero values of Wi are determined. 

As an example, in the forward case the singular values for all GFFs other than A1^q{0) are zero and the result is 
the minimization in the subspace of A'q(O) alone. In this case, A^g(O) is still determined simultaneously from several 
different lattice measurements and this procedure thereby automatically incorporates more physical constraints than 
the conventional method of fitting a single lattice observable. 

The errors in J-' are calculated using a jackknife analysis in the following way. The system (jl8|l is solved simul- 
taneously for the full set of lattice measurements R to get the average solution vector J^. Then, the appropriate 
subsamples Ri are formed from which the subsample solutions !Fi are calculated. The error vector for the solution, 
CTjr, is formed from the J-^iS by the familiar jackknife formula (see e.g. |l8l | and references therein for details). 

In summary, the method introduced in this section provides a way to extract maximal information from a set of 
lattice calculations. It simultaneously includes all available data in a unified manner and corrects for correlations 
present in the underlying sample of gauge field configurations. 



III. RESULTS 



In this section, we apply our procedure to the case n = 2 using the operators in Eqs. (jA6IA7|l and thoroughly 
analyze the results. For this exploratory calculation, we use SESAM unimproved Wilson configurations on a 16'^ x 32 
lattice at P = 5.6 with Kgca = ^vai = 0.1560. We note that this is the heaviest of the three quark masses used in earlier 
calculations of moments of quark distributions [T3| . and using the lattice spacing = 2.01 GeV from the chirally 
extrapolated nucleon mass, the pion mass (am),r = 0.446(3) corresponds to niT^ = 896(6) MeV. All the computational 
details for calculating smeared sources and sinks, nonrelativistic sequential sources, the difference approximations to 
the operators C'^^^^^ ^ }' ^P^'^ polarization projection, and Dirichlet boundary conditions are as in Ref. [l^ . In this 
present work, we relabel the time slices such that the source is at Tgrc = 1, the sink is at Tgnk — 13 and boundaries 
are at t = —9 and t = 22. 

The one-loop perturbative renormalization constants for the cases we consider in this work are taken from Ref. |14| . 
Note, that in the case of the operators O^^^^, Eqs. (|A6IA7(I . the renormalization constant will depend on the index 
combination chosen. The diagonal and non-diagonal index combinations belong to distinct irreducible representations 
of the lattice point group and have different renormalization factors: ^diag(A'^ = 4GeV^) — 0.9768 and ■^non?diag(M = 
4GeV^) = 0.9884. 



A. Comparison of Minimally Determined and Overdetermined Measurements 

To explore how the method works and to compare it with the conventional method of measuring a single operator 
for each observable, we consider the case of virtuality iviow — — 0.5925 GeV^, which may be achieved by choosing any 
of the momentum combinations displayed in Tab. ^ 

We start by considering only the combination of external momenta, P ' = (0, 0, 0), P = (— p/, 0, 0), listed in the first 
column in Tab. The system of equations, Eq. H15|l then contains the operators 

(J^'|OS[ag,ll^> ■ (^PS[ag,2l^> , (^'|OSifg,3l^> . (^Puo}IO , (^'|0{20}l^) - (^'|0{21}l^) , (20) 



with the operator index conventions given in Appendix El Eqs. HA6IA7|) . 
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TABLE I: Different combinations of external momenta with virtuality fviow = — 0.5925GeV^. 



N 


1 


2 


3 


4 


5 


6 


7 


P' 


(0,0,0) 


(0,0,0) 


(0,0,0) 


(0,0,0) 


(0,0,0) 


(0,0,0) 


(-Pi, 0,0) 


P 


(-Pi, 0,0) 


(0,-Pi,0) 


(0,0, -Pi) 


(Pi,0,0) 


(0,P!,0) 


(0,0, 


(0,0,0) 



Plateau plots of the ratios, R{t, P' , P), are shown in Fig.Hlfor three matrix elements, (-P'|CSiag,il-P> ' (-f"Pdiag,2l-P) > 
and (P'|C'|2o}l-P')- Measurements are averaged over the interval [rmin = 5,rmax = 9] and this interval is denoted by 
vertical lines in the plots. The solid horizontal lines show the extracted plateau values and the dotted lines denote 
one standard deviation. 

The resulting lattice matrix elements, together with their jackknife errors are as follows: 

^d4i(^''^) 0.2327 ±0.0224, 

^d4.2(-f"'^) = 0.2970 ±0.0235, 

i?d;fg_3(P',P) = 0.00832 ±0.01651, 

i?{;o}(P',P) = -0.09768 ±0.01648, 

^{20}(-f"'^) -0.1802 ±0.0205, 

^{21} (-P',^) = -0.05926 ±0.01162. (21) 

We now consider a minimal set of three operators to determine, but not overdetermine, the three generalized form 
factors, ^25"^, 520"^, and 6*20'^ • The subset of diagonal operator equations alone are insufficient T^, since the diagonal 



^20 ' ^20 ' '^'-^'^ ^20 

index combinations are linearly dependent and provide only a system of equations of rank two. 

The minimal choice to determine, but not overdetermine, the form factors is to choose two diagonal operators and 
one off-diagonal operator, analogous to familiar procedures to calculate moments of parton distributions 10, 14] or 
electromagnetic form factors |l3l 12^. The best choice for this purpose is the set of matrix elements (P'lOdjag il^*), 
(-P'l^diag,2l-P)' and (P' 1 0^20}!^)' "^^ich have the smallest relative errors and provide a system of three linearly 
independent equations. 

The form factors and associated errors determined this way are shown in Fig.|21by the triangles plotted to at the 
left at = 0. 

Still considering the same single choice of external momenta, the full set of six matrix elements in Eq. (j20|l provides 
an overdetermined set of equations that yield the form factors and errors labeled = 1 in Fig |21 We previously 
argued that the C form factor is least accurately determined because its coefficients have the highest powers of A, 
and the minimally determined fit, N — bears out this expectation. Note that the overdetermined fit, 1 leaves 
the errors in A and B roughly the same, but substantially reduces that of C. The results of the two fits agree within 
errors, but since the overdetermined fit exploits the maximal information in the lattice measurements, Eq. H21|) and 
has smaller overall statistical errors, we consider it superior. 

We now progress to the next level by including all the external momentum combinations enumerated in Tab.Q] The 
results are displayed in Fig. [21 where the abscissa denotes the total number of momentum combinations included in 
the fit. Entry N corresponds to the inclusion of momentum sets 1, 2, • • • of Tab.Q] As expected, the errors for each 
form factor decrease significantly as new statistical information is included by adding momenta in new directions. 
The improvement is generally less than ^/N because measurements of different momenta on the same lattices are 
correlated. 

The overall result in Fig. |21is a strong validation of our approach. Comparing the minimally determined, A^ = 
result at the far left with the overdetermined result with all seven momentum combinations, A^ = 7 at the far right, 
we observe that our method reduces the error in the least accurately determined form factor C by a factor of 5, the 
error in B by a factor of 3, and the error in A by a factor of 2. 



B. Generalized Form Factors and Quark Angular Momentum 



Having validated our analysis technique, we now apply it to the full set of external momenta and virtualities listed in 
Appendix^ Tab.^l The results of the full, overdetermined analysis for the flavor non-singlet generalized form factors 
^20"^ (0; ^20^ (i) J and C2~'^^{t) are shown in Fig.O We observe that with the present analysis, all three generalized form 
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FIG. 1: 



Plateau plots of the ratios R{t, P' , P) for the n = 2 operators Odiag.i, C'di'^ 



ag,2i and C{20}- 



factors are determined quite accurately throughout the full range of virtuality up to <: = 3.1 GeV^. As expected, some 
combinations of external momenta induce more statistical noise than others, but the overall structure of the form 
factors is well determined. Although there is no fundamental argument for the functional dependence on i, it is useful 
to fit the form factors to the dipole form that provides a good phenomenological fit to the nucleon electromagnetic 
form factors, and the results of least-squares fit by dipole form factors are shown for A{t) and B{t). Although the 

individual up and down form factors C2^'^{t) are non- vanishing, their difference, C2~'^{t), is statistically consistent 
with zero. We note this flavor independence is a feature of the Chiral Quark Soliton Model in Ref. j2ll |. 

Although results are not yet available at other quark masses to enable extrapolation to the physical quark mass, 
these results clearly show the behavior of the generalized form factors in the 900 MeV pion world. For this non-singlet 
case, there are no corrections from disconnected diagrams. 

It is useful to point out that in all cases for which i ^ 0, there were no zero (or nearly zero) singular values 
in our calculations. Hence, in these cases, the singular value decomposition is completely equivalent to minimizing 
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FIG. 2: Generalized form factors obtained by simultaneous fits to A'^ external momentum combinations having virtuality tviow 
As described in the text, the A*' = points, denoted by triangles, use three operators at a single external momentum combination 
to determine the three form factors. The remaining points, denoted by squares, use six operators and A'^ external momentum 
combinations to determine the three form factors. 



Eq. (|18|l by conventional least-squares analysis. In the case of i = 0, because of the explicit factors of A appearing 
in Eq. ((SJ, only the coefficients A^jq(O) can be determined. As expected, the singular value decomposition handled 
this automatically, with the zero singular values resulting in minimization in the appropriate subspace of A'q(O). 
Note that this procedure introduces no bias or model assumptions. Rather, it simply specifies the correct physical 
subspace. 

Since the total quark angular momentum is given by the zero virtuality limit, Jq — ^[^20^'^(0) + i?2o^'^(0)], it is 
particularly interesting to study this limit. Two issues need to be addressed. The first is extrapolation to t = 0. This 
is not a problem for A, since it is just the moment of the spin-averaged parton distribution j42o(0) = {x). However, 
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FIG. 3: Generalized form factors A2o^{t), 820^ (t), and 6*2''^ (t) for all available virtualities obtained using the full set of operators 
and external momentum combinations. The dashed curves denote dipole fits to A and B to guide the eye and extrapolate 
to t = 0. The form factor C is consistent with zero. Four data points, denoted by horizontal brackets, have been shifted by 
\t\ — 0.2 GeV^ to the right for clarity in plotting. 



like the electromagnetic form factor F2, B cannot be measured at i = since the kinematical factors always contain A, 
and therefore it must be extrapolated. From the behavior of the flavor non-singlet combination i?2o'^(t) in Fig.|3| one 
might expect significant extrapolation errors. However, as shown in Fig. 01 the cancellation between the up and down 
quark contributions is nearly complete and the flavor singlet combination -620^'^ (i) is consistent with zero. Hence, 
essentially the full contribution comes from A''2q'^{Q). 

The second issue is the fact that disconnected diagrams, which are beyond the scope of the present work, must also 
be included in these flavor singlet matrix elements. 

It is interesting to note that light cone arguments indicate that the contribution of each Fock space sector to 
the gravitomagnetic moment B{Q) — Yl, q g ^213 i*^) = is separately zero Although this may not necessarily 

imply that connected and disconnected diagrams for i?(0) separately cancel, it is suggestive of our result above that 
the connected contribution to i?(0) was consistent with zero. Note that the separate vanishing of quark and gluon 
contributions to B{Q) was also conjectured in Ref. p^ . 

With the caveat that we must omit disconnected diagrams for the present, we proceed to quote our lattice results 
from connected diagrams. Using the results of Ref. 1^1 at k = 0.1560, the quark spin contribution to nucleon angular 
momentum is given by the lowest moment of the longitudinal spin distribution 

^AE - \[{l)Au + {l)Ad] 

~ i 0.682(18). (22) 
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bmce Bw'^iO) ~ 0, the total quark spin and orbital contribution to the nucleon angular momentum is 



Jn = 



B. 



u+d 



1 

2 

- 0.675(7). 

2 ^ ^ 



(0)] 



(23) 



Thus, since these two results are equal within errors, if we consider only connected diagrams in the world of 900 MeV 
pions, we conclude that the quark spin produces 68 percent of the nucleon angular momentum and the quark orbital 
angular momentum contribution is negligible. Although this is superficially consistent with expectations from rela- 
tivistic quark models, we note that serious physical interpretation and quantitative comparison with quark models 
requires consideration of the disconnected diagrams. The striking fact for our present purposes, however, is that it is 
possible to measure the connected part to the order of one percent, under the assumption that B^jn'^jO) = ± 0. 

Note that quenched connected diagram contributions of AE and Jq have been reported in Refs. (23. l25l|. Since 
those calculations were performed for pion masses ranging above 500 MeV, we would expect quenching effects to be 
negligible and hence that their results would be statistically consistent with ours. The chiral extrapolated results of 
Gockeler et. al. that 66 ± 14 percent of the nucleon spin arises from the quark spin and that 6 ± 14 percent arises 
from quark orbital angular momentum are completely com pat ible with our results because of the small quark mass 
dependence. The results for AS by Mathur et. al. in Ref. |2J| are also consistent with ours. Reference [lj| showed 
that the linear extrapolation of AE is essentially constant, only changing from 68 ± 2 percent of the nucleon spin at 
TOtt = 900 MeV to 69 percent in the chiral limit. This is consistent with the connected diagram chiral extrapolation 
by Mathur et. al. of 62 ± 8 percent. The results for Jq, however, are inconsistent. From the results of their Fig. 3 
and their 1.045 renormalization factor, Jq contributes 84 ± 8 percent of the spin, which disagrees with our result of 
68 ± 1 percent. We believe the discrepancy arises because Mathur et. al. use the dipole prescription to extrapolate 
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the sum A2q'^ + -820^"^, whereas, as shown in our Fig.01 we calculate A2q'^{0) directly without extrapolation and only 



^20 ^ ^20 

-"20 



extrapolate the nearly vanishing ^2^"^ 



IV. CONCLUSIONS 

In summary, we have presented a new method for calculating generalized form factors in lattice QCD that extracts 
the full information content from a given lattice configuration by measuring an overdetermined set of lattice observ- 
ables. We demonstrated its effectiveness in an exploratory calculation of 71 = 2 generalized form factors up to 3 GeV^ 
and showed that it reduces errors to as small as one-fifth of those obtained by a conventional minimally-determined 
analysis. The final error bars for A2^(t), A^^'^it), and -620'^ (i) are typically of the order of five to ten percent, provid- 
ing useful information about the Fourier transform of the transverse structure of the nucleon. Because of fortunate 
cancellation in the fiavor singlet case, the connected diagram contributions to the total quark angular momentum are 
measured to the order of one percent. Although this exploratory calculation was performed for the heaviest of the 
three SESAM quark masses used in earlier calculation of moments of parton distributions, we expect the technique 
to be sufficiently robust to treat all the masses. 

We note that the first calculation of generalized form factors were reported by Schrocrs in 19] using the standard 
minimally-determined analysis, and more extensive calculations by the QCDSF collaboration using our method are 
being published simultaneously with this present work [25l | . 

This work provides the foundation for a number of promising investigations. The n = 3 relations introduced in 
this work enable us to calculate the i-dependence for three moments and thus explore the variation of the transverse 
structure of the proton with x. The same methodology can, of course, be used to calculate spin-dependent generalized 
form factors. In the longer term, extensive calculations with emerging multi-Teraflops computers dedicated to lattice 
QCD will enable the chiral and continuum calculations required to have definitive impact on the new generation of 
experiments currently being undertaken . 
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APPENDIX A: MATRIX ELEMENTS IN TERMS OF GENERALIZED FORM FACTORS 

In order to find the relation between Minkowski-expressions, e.g. ^^m^ - ^j^^ corresponding Euclidean 

quantities, (O^)^^ ^ ^ j, we observe that the operators under consideration are constructed from gamma matrices 
and covariant derivatives D. The conventions we use, consistent with Ref. are 

4 

h""), = E '^M^ (^"^y ' '^M^ = diag (i, i, i, 1) , 
i/=i 

4 

i^"")^ = EV(^'')^ V = diag(l,l,l,-i) , 

and we will set (a^)"-" = (a^)''-^ Since all operators are symmetrized, the general transformation rule for the 
operators with one gamma matrix and n — 1 covariant derivatives investigated in this work is given by 

(e^") - h,,^A.^, . . . (o-)^''--''"^ , (Al) 

where a summation over the indices i^i is implicit. In the following, we will drop the label E or M, and it will be clear 
from the context if Euclidean or Minkowski expressions are used. 
We write matrix elements in terms of Dirac spinors as follows 



{P'\0\P) =U{P')JCo{P\P)U{P). 
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The ratio P) can be written explicitly in Minkowski space as 

I ground state 



^a"'(^',^)|,.o.nd.tatc - {E{P')E{P)y''' 



1 1 ^ -^'^ 

(/P' + m)]-Tr[r„„poi(/P + m)] 

xiTr[rpoi(/P' + m)/Co(P',P)(/P + m)] , (A2) 

neglecting all excited states which would introduce a time dependence. Note that i?Q"*(P', P) is directly proportional 
to {P'\0\P). The ratio in Eq. ijSSIl is expressed in a general form that applies to all quark bilinear operators and 
projectors Fpoi/unpoi relevant to this work. 

The Minkowski space expressions for the QCDSF and the LHPC projectors used in our analysis are as follows 

^\n^S = ^ (1 + 70) , r^HPC ^1(1 + (1 _ ^^^3) ^ (A3) 

^T^r - ^ (1 + 70) , r^^-- = 1(1 + ,o) 7572 . (A4) 

We also record the complete expressions for the kernels ICo{P' , P) in Eq. (|A2|) . Following the general parameterization 
of 0, we have for the lowest three moments rt = 1, 2, 3 in Minkowski space 

(p'l^iMlA's...} |p^) ^ UiP')JCl^'"''-\P',P)U{P) 

{P'\i^^r'i;,\P) = I7(F')|7^^A.i(A^)+i^^^S,.i(A')}c^(^) 

A{A'iAA'2} ^ 

+ H) C,,2(A2) C/(P) 



{P'\ij^j^^^'D^^'D''^Hg\P) = U{P')[{~i)'P^'''P'"'r^^A,.so{^') 

+ (-i)2A{^iA^=7A'3}A,32(A2) 



i(-i)2A{^^A^=^^-^B,,32(A2)| C/(P) . 



(A5) 



In practice, we evaluate the ratio - ^ (p' by inserting the kernels JC^^^^'"^ from Eq. l|A5p and the 

projectors in Eq. (|A3IA4|I into Eq. I|A2|I . For given values of the nucleon momenta P and P', the ratio is computed 
numerically, determining the kinematical coefRcients of the GFFs A{t), B{t) and C{t). Eventually, we transform the 
free indices {/ii/i2 . . . } to Euclidean space with the use of Eq. (|A1|) . The system of linear equations is then constructed 
by equating the Euclidean lattice result and the transformed Minkowskian parameterization. 



1. Lattice operators 

On the lattice side, we have H(4)-invariance instead of Lorentz-invariance, which complicates the classification and 
the renormalization of the operators. Following 1301] . we choose appropriate linear combinations of the operators 
^ tj.2 ■ ■ ■ } ' corresponding to representations of H{A). 

For the case of one derivative, n — 2, there are three diagonal operators 

of "^''^^^ = \ [On + O22 - O33 - O44] , 



o^''^-'^^' = ^[^33-04 



4''''''^' = ^[011-022], (A6) 
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which are symmetric and traceless by construction. AdditionaUy there are six non-diagonal combinations 

0no„^diag,„=2 ^ 1 [(^^^^^^ ^ 0^^^^] ^ 0^^^^^^^ = 1 . . .4; /.I < /i2 , (A7) 

giving in the most general case a set of nine independent operators. There is no operator mixing for n — 2. 
In the case of two derivatives, n = 3, we take 

n— 3 I I [/n /o 1 /nil— 3 



or' = ( I) [0{2ii} - 0{233}] , or' = ^ 



[C{122} - 


(- C'{133} 


- 20{i44} 


[C{211} - 


1- C'{233} 


- 20(244} 


[C{311} - 


1- 0{322} 


- 20{344} 


[C{411} - 


1- C'{422} 


~ 20{433} 



or' ^ [^{311}-0{322}], 06='-2l/2 

or' = (^) ^ [Oi,n} - 0{422}] , ^8=' = ^ [^^{411} + 0{422} ^ 20{433}] , (A8) 

while the four non-diagonal combinations arc denoted by 

/-/-)non-diag,n=3 ,, i A. ,, ^ ,, ^ ,, 

Ot,u^^2A = <-^{miA'2A'3}' /il,2,3 = 1 • ■ -4, Ml < Ai2 < Ai3 , 

giving all in all 12 traceless and symmetric linear combinations. It turns out that some of the operators in Eqns. (|A8|I 
mix with different non-symmetric representations (cf. Ref. jHSl) under renormalization. Fortunately, the mixing 
coefficient turns out to be negligible, at least for the specific case considered in iXQj. It has to be checked if this holds 
true for all potential mixing candidates in the n = 3 case. 

The choice of index combinations for the construction of the H(4) operators is not unique. However, in the continuum 
limit all options will lead to the same GFFs. The interesting question of determining an optimal set of operators that 
minimizes lattice artifacts is beyond the scope of this presentation. 



2. Lattice momenta 



The possible values of the nucleon momenta P and P' and therefore the momentum transfer squared t are severely 
restricted by the lattice momenta. The general three-momentum on a periodic lattice is given by 

— > 2n 2tt 

P = = ~~r (Pi'P^'Ps) ' P^ = -Ls,-Ls + l,...,Ls-l,Ls, 

aLg aLs 

where a is the lattice spacing in GeV~^ (we set a = 1 in all intermediate steps), and Lg is the spatial lattice extension, 
which is in our case 16. It is clear that the lattice spacing times the spatial extension determines the smallest possible 
non-zero momentum, while the lattice spacing alone fixes the largest available momentum. In practice, however, we 
are even more restricted to values pi <^ Lg, because large momenta lead to considerable noise. For the nucleon four 

momenta P and P' , we use the continuum dispersion relation, Pq = P4 = \J m"^ + P'^, where m is the nucleon mass, 
determined in our lattice calculation. In order to construct the maximally determined set of linear equations, Eq. H14|l . 
we need to know all leading to the same momentum transfer squared t. So far we have propagators available 

for the two sink momenta — (0, 0, 0) and = (—1, 0, 0), and the value for ~q = ~p' — is restricted by < 6. 
We plan to extend the allowed values in future investigations. The total number of different t is then 16, including 
the forward case. 

Note, that t depends in a nontrivial manner on and ^ individually through the energies i-4,P4, and not only 
on the relative momentum ~q . For 8 of the 16 possible values of t, we encountered highly visible fluctuations in the 
nucleon propagators with negative values for the two point functions at large t ^ Tgrc- We excluded the corresponding 
momentum combinations from our analysis. A list of all remaining ~q^,~p*',~p' and t, is shown in table Hll 
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